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Abstract. We investigate the statistical equilibrium properties of a system of classical particles interacting 
via Newtonian gravity, enclosed in a three-dimensional spherical volume. Within a mean-field approxi- 
mation, we derive an equation for the density profiles maximizing the microcanonical entropy and solve 
it numerically. At low angular momenta, i.e. for a slowly rotating system, the well-known gravitational 
collapse "transition" is recovered. At higher angular momenta, instead, rotational symmetry can sponta- 
neously break down giving rise to more complex equilibrium configurations, such as double-clusters ("dou- 
ble stars"). We analyze the thermodynamics of the system and the stability of the different equilibrium 
configurations against rotational symmetry breaking, and provide the global phase diagram. 

PACS. 05.20.-y Classical statistical mechanics - 04.40.-b Self-gravitating systems - 64.60.Cn Order- 
disorder transformations; statistical mechanics of model systems 



1 Introduction 

In this article, we study the equilibrium properties of a 
system of N classical particles subject to mutual gravita- 
tion, assuming that this self-gravitating gas is enclosed in a 
finite three-dimensional spherical box and rotates around 
its center. The Hamiltonian reads 

1 ^ 

EE H^an}, {pj) = —J2p^+ <Z.({r J) (1) 



(2) 



l<i<j<N 



where € C M'^, G M'^ and m > denote, re- 
spectively, the position, momentum and mass of the i-th 
particle, while G is the gravitational constant. (In the fol- 
lowing, we set m = 1.) V stands for the (volume of the) 
box containing the particles. Notice that the total poten- 
tial energy scales as iV^. One wants 

a. to calculate the spatial distribution of particles at equi- 
librium, i.e. the most probable microscopic state; 

b. to derive the global phase diagram of the system; 

c. to study its thermodynamics (e.g. caloric curves); 

d. to describe the phase transitions that eventually take 
place. 

In general, the physics of systems whose microscopic 
constituents interact via long-range potentialg|^ such as 

^ By which we will mean decaying with the interparticle dis- 
tance r more slowly than r'^'" with e > in _D dimensions 
when r oo. 



(|l|,y|), is highly non-standard and the analysis of their equi- 
librium state represents a considerable theoretical chal- 
lenge At odds with more conventional systems with 
short-ranged forces, they are not additive (i.e. they can- 
not be divided into macroscopic subsystems with negligi- 
ble mutual interaction) and not extensive (i.e. the den- 
sities of thermodynamic functionals are not bounded in 
the limit N — > oo). These facts have several important 
consequences. Long-range systems can have negative spe- 
cific heat (^-|^; statistical ensembles can be inequivalent 
[^[t HioIFI ; they do not possess a proper infinite- volume 
limit P 1[ ; and they can attain inhomogeneous configura- 
tions (indeed, even their ground state is inhomogeneous, 
and it has been suggested that fractal structures may 
as well arise ^). Conventional statistical mechanics 
techniques that apply to homogeneous, short-range sys- 
tems hence fail for long-range systems. The key theoret- 
ical problem here is a very fundamental one: to devise a 
mathematical framework that allows the study of phase 
transitions and other collective phenomena . 

Self-gravitating gases are the most prominent exam- 
ple of a long-range system and have attracted a great 
deal of attention over the years. Most static theories rely 
on the microcanonical ensemble, where conserved macro- 
scopic observables represent the natural control param- 
eters (see for a review). The task in this setting is 
that of finding the most probable (entropy-maximizing) 
equilibrium configuration of (|l|,H) with V finite as a func- 



In fact, the paper |g] pointed explicitly to the failure of the 
canonical ensemble near first-order phase transitions in gen- 
eral, and to its non-equivalence with the fundamental micro- 
canonical ensemble, which displays a negative heat capacity 
there. 
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tion of the integrals of (Hamiltonian) motion, the simplest 
and physically most relevant being the total energy E and 
the total angular momentum L — X^i^i '"i ^ P^^ whose 
conservation is related to the invariance of (|l|) under ro- 
tations. Taking the conservation of L into account leads 
however to enormous technical difficulties. Most authors 
have thus neglected rotation by breaking the rotational 
symmetry explicitly from the outset, e.g. by taking a non- 
spherical V, and used E as the only control parameter 
(see e.g. |16| - [l8| for some recent work). The analysis of 
the equilibrium state of self-gravitating gases is hence a 
rather well-studied problem without rotation. It turns out 
that at high energy (i.e. temperature), where the kinetic 
term dominates, the system is most likely found in a ho- 
mogeneous cloud (shortly, "gas" state) filling the available 
volume. At low energy, instead, where the gravitational 
energy dominates, a collapsed configuration is preferred, 
where particles form a single dense globular cluster lying 
in an almost void background ("single star"). This is the 
well-known gravitational collapse transition first described 
in In between these two "phases" (not being homo- 
geneous, the single-cluster is not a proper thermodynamic 
phase), for a whole range of energies, the specific heat is 
negative. In this transition regime the canonical ensemble 
fails|. 

For a rotating system, the situation is expected to 
be substantially more complex. From a qualitative view- 
point, the equilibrium density profile will depend on the 
ratio between the rotational and gravitational contribu- 
tions to the total energy. When the latter dominates, grav- 
itational attraction should cause the system to collapse. At 
high ratios, instead, when rotation is sufficiently fast (high 
angular momentum), more complex distributions should 
arise. Dynamical studies based on fluid-mechanics tech- 
niques [pi|-23 suggest that ring- like and disk-like struc- 
tures might appear. Ultimately, at sufficiently high rota- 
tional energies, two distinct dense clusters (i.e. a "double 
star") should form. 

The richness suggested by the fluid-dynamical picture 
cannot be recovered in a static equilibrium theory with- 
out the inclusion of rotation. Double-cluster configura- 
tions can arise in a static framework only from the spon- 
taneous breaking of the rotational symmetry of (|l|) , which 
should take place when the angular momentum is suffi- 
ciently high. Effects connected to rotation should also lead 
to the formation of rings and other types of structures. 
Despite some attempts P^,p5|, however, a detailed static 
theory embodying angular momentum is lacking. 

The purpose of this work, which builds on ^ , is 
to include angular momentum in the microcanonical the- 
ory. Using a mean-field approximation, we derive an inte- 
gral equation for the density profiles corresponding to sta- 
tionary points of the microcanonical entropy surface and 
solve it numerically as a function of E and L = \L\. The 
usual collapse transition is recovered at low angular mo- 
mentum. At high angular momenta, instead, we find that 



a spontaneous breakdown of the rotational symmetry oc- 
curs at sufficiently low energies. This gives rise to more 
complex equilibrium structures, including "double stars", 
rings and disks. We derive the global phase diagram of the 
self-gravitating gas in the {E, L) plane. By studying the 
Hessian of the microcanonical entropy ^, we charac- 
terize three pure phases (where the system is found in a 
"gas", "single star" and "double star" configuration, re- 
spectively) and a large mixed phase with negative specific 
heat, phase separation and competition between differ- 
ent equilibrium density profiles. Finally, we analyze the 
thermodynamics of the system, deriving the caloric curves 
(temperature versus energy) in the different phases, and 
analyze the stability of the stationary distributions. 

This paper, which follows |^^, is structured as follows. 
In Sec. 2 we expose the microcanonical mean field theory 
of (H), derive the entropy functional and the stationar- 
ity condition. Sec. 3 is dedicated to the results. We report 
the numerical solution of the main equation, together with 
the global phase diagram and a few equilibrium configu- 
rations. Then we pass to the thermodynamics of the sys- 
tem, with special emphasis on the rotational-symmetry- 
breaking transition, the physics of mixed phase, and the 
stability problem. Finally, Sec. 4 contains our conclusions 
and a some remarks about the work presented here, and 
a list of open problems. 



2 Microcanonical mean-field theory 

We consider the system with Hamiltonian as given in (|l|,| 
enclosed in a three-dimensional spherical volume V , to 
preserve rotational symmetry and ensure conservation of 
the total angular momentum L. At the same time, the 
box breaks translational invariance, hence the total linear 
momentum is not conserved. The aim of the microcanon- 
ical theory is to find the particles' density profiles p{r) 
satisfying 



p(r)dr — N 



(3) 



It has been recently shown that ensemble equivalence is 
restored in the "dilute" limit (iV, V) cxi with TV/l/^/^ fixed, 
in which thermodynamic functionals exist |20|. 



that maximize the entropy {k = 1) 

SNiE,L) ^liiWNiE,L) (4) 
Wn being the microcanonical "partition sum" (/i = 1) 

Wn{E,L) = ^ / S{Hn - E) 5{L -J2r,x p,) Dr Dp 
■ i=i 

(5) 

We used the shorthand notations Dr = Hili ^'"i 

= Yli=i'^Pi- e is a constant that makes Wn dimen- 
sionless. Integrals over momenta are from —oo to -|-oo, 
while integrals over {r^} are performed over . Clearly, 
such equilibrium profiles will depend on E and L. 

We now calculate the microcanonical partition sum 
(^. To perform the integrals over momenta, namely to 



E.V. Votyakov, A. De Martino, D.H.E. Gross: Thermodynamics of rotating self-gravitating systems 



3 



evaluate 

FN{{r^},K,L) 



N N 

= / -^Y.P'^ '^(i: - ^ n X Dp (6) 

one can follow Laliena |^ (see also |^ ) and use the Laplace 
transform of Fn in K, that is (5fts > 0) 



Fjv({r J, s,L)^ / FN{{r,}, K, L) e 
Jo 



-'^ dK 



(7) 



Inserting the integral representation of the 5 function and 
performing the trivial integral over K one gets 

FN{{ri\,s,L) ^ j e''^--^-'5:,'^-r,xp,-f E.p? 

(8) 

where Doj = da;/(27r)'^. The above integrals are at most 
of Gaussian type and can be performed to yield 



FAr({r,},s,L) = 



/n N 3JV-3 i„r^i-ir 

(27r)— 2— e~2^-^ ^ 



Vdetl 



(9) 



S 2 



where I = I({ri}) denotes the inertia tensor, with ele- 
ments (a, 6=1,2, 3) 



N 



Iab{{ri}) = ^{rfSab - ri^an,b) 



(10) 



and L'^I = Lal^^Li,. The inverse Laplace trans- 
form of (0) is given by 



FNi{n},K,L) 



1 



^ T 1 , 3JV-5 

r(^)Vra(^-2^^^^)^ 

(11) 



for K > \L^l-'^L and Fjv({r J, if, L) = otherwise. 
Hence after integrating out the momenta the microcanon- 
ical partition function reads 



N\J Vd^ 



Dr 

(12) 



where A = (27r)^— /r((3iV - 3)/2). The term in square 
brackets in (O) is nothing but the kinetic energy of the 
system. Now setting for simplicity JC — E — ^lFi^^L — 
<?({ri}), we remark that the integrand is 



gAf[f log/C- 2^ log/C-jJ^ log \/d5tl 



(13) 



Being interested in the behaviour of the system for large 
N (see below), we shall retain just the leading order in N 
in the above expression. We are thus left with 

W^^(£;,L) = ^ j[E-\L'^l-^L-1>{{r.,})]'^ Dr 

(14) 



and it remains to integrate over . 

To this aim, we write the potential and the com- 
ponents of the inertia tensor as functionals of the density 
profile p as follows: 

<^({r-J) - m = -^J ^yzVl-drdr' (15) 

/ah({r.}) ^ Iab[p]= I p{r){r^5ab~ran)dr (16) 



Notice that in this way two- and many- body correlations 
are neglected. This allows to recast ( [l7| ) in the form of the 
functional-integral 

(17) 

(mf = mean field) where P[p\ is the probability to observe 
a density profile p = p(r). To estimate the latter, we follow 
the logic of Lynden-Bell . We subdivide the spherical 
volume V into K identical cells labeled by the positions of 
their centers. The idea is to replace the integral over 
with a sum over the cells. In order to avoid configurations 
with high densities where other physical processes (e.g. 
nuclear reactions) become more important than gravity, 
and cure the short-distance singularity of the Newtonian 
potential, we assume that each cell may host up to riQ 
particles (1 <C no ^ N). This condition is essentially 
equivalent to considering hard spheres instead of point 
particles. P[p\ is now proportional to the number of ways 
in which our TV particles can be distributed inside the K 
cells with maximal capacity no. Denoting by n(rfc) the 
number of particles located inside the /c-th cell, a simple 
combinatorial reasoning leads to 



P[P\ 



n(ri 

m n 



n 



no! 



no 



cells k 



. \n{rk) 



(18) 



where it is understood that the product involves configu- 
rations {n(rfc)} such that n(rfc) = N. Introducing the 
relative cell occupancy 



c(r) 



n(r) Vp{r) 
Kno 



no 



(19) 



and approximating the factorials by means of Stirling's 
formula (assuming n{rk) ^ 1 and no — n{rk) ^ 1), we 
get 

P[c] CX A^! Q-^^ I [^(r)logc{r) + il-c{r))log{l-c{r))]dr ^ 

= A^! e^'§'-'^t'^(^)'°s''('^^+(^^^("'^^'°s^^^^^'^^)l'''^ (20) 

where we introduced the dimensionless variable x — r/R 
and defined the average coverage 







NV 
noKR^ 



c{x)dx 



(21) 
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It is simple to check that in terms of c{x) the potential 
and the inertia tensor are respectively given by 



c{x)c{x') 



dxdx' 



Iab[c] = 



c{x){x Sab ~ XaXb)dx 



(22) 



(23) 



In the following, we shall measure energies in units of 
and inertia tensor components in units of NRl^ so that we 
will be dealing with the reduced (dimensionless) quantities 



•PI 



1 

202 



c{x)c{x') 
\x-x'\ 

„2) 



dxdx' 



(24) 



Iab[c] = J C{x){x 6ab " XaXb)dx (25) 

Plugging (^0|) into (p7|), one arrives at the familiar form 
W^He, L)o^ j e^'^" W dc{x) (26) 
where the "action" S^f has the following expression: 



S 



mf r 

N V 



log 



Ttt-I 



L - <?[c] 



y"[c(ic) log c(a;) + (l-c(a;))log(l-c(a;))] dx 

(27) 

Clearly, I = I[c]. S^f has an obvious physical interpreta- 
tion as the sum of the energetic and combinatorial con- 
tribution to the entropy, respectively. For large N, the 
integral ( ^ ) can be computed as usual by the steepest- 
descent (Laplace) method. This immediately yields 



W^\E, L) ~ exp 



7VmaxS']^f[c] 

c(x) 



(28) 



hence the "physical" value of the entropy density for large 
N is nothing but the maximum of S*^^ over the space of 
relative cell occupancies c. 

An elementary variation of with respect to c, the 
constraint on being enforced by a Lagrange multiplier 
playing the role of a chemical potential, straightforwardly 
leads to the stationarity condition 



log- 



c{x) 



iuix) 



X xf 



/i 



(29) 



1 - c{x) 
or, equivalently, 

where u; = w [c] is the angular velocity (related to the total 
angular momentum by the relation i = lu)), and /3 = /3[c] 
and U{x) are respectively defined as 



(30) 



/3 = 



3/2 



[E ^ 



c{x') 



U{x) = - I ^^(^ dx' = 202 
J \x-x'\ 



2 S'P 



5c{x) 



(32) 



One sees that (3 is related to the kinetic energy of the 
systems, i.e. to the (inverse) temperature. The essence of 
the mean-field approximation is clearly expressed by the 
fact that 



(33) 



Equatio n (29[ ) (or (|3C|)) is our central result. Functions c* 
solving (E9) and being entropy maxima in the space of c's 
represent our desired equilibrium distribution of particles. 

The correct way to analyze the problem consists in 
solving ( p^ ) at fixed energy and angular momentum, sub- 
sequently calculating intensive quantities (temperature and 
angular velocity). For the sake of simplicity and without 
any loss of generality, we now fix the angular momentum 
to lie parallel to the 3-axis, and concentrate on \L\ = L. 
We remark at this point that Lynden-Bell statistics, which 
is reminiscent of Fermi-Dirac statistics in real space (i.e. 
not in phase space), plays a crucial role. In fact, once over- 
lapping is ruled out, the Hamiltonian (|l|) has a well-defined 
ground state, with particles collapsed in a core but with- 
out coming too close to each other. If one used Boltzmann 
statistics and point particles, instead, the system would 
have no ground state, since the potential energy would be 
unbounded from below. Antonov catastrophe p9| | can be 
seen as a direct consequence of this fact. For this reason, 
Lynden-Bell statistics is probably more appropriate for 
self-gravitating systems^. The existence of a ground state 
ensures that ( p9| ) will always have a solution at fixed E 
and L. Of course, there may be multiple viable solutions 
at the same E and L, each bearing its entropy. In such a 
case, the criterion is simply that the higher the entropy, 
the more probable the solution. 

Upon varying E and L, one can explore different re- 
gions of the parameter space and ultimately obtain the 
global phase diagram in the whole {E, L) plane. The main 
effect one expects from the inclusion of rotation is that, 
for sufficiently high angular momenta, upon decreasing the 
energy from high values, rotationally-symmetric solutions 
(e.g. homogeneous clouds) will become unstable against 
fluctuations that break rotational symmetry, and solutions 
without rotational symmetry (e.g. "double stars" ) will bi- 
furcate continuously from them. 

The main problem at th is p oint is merely technical. 
One can only hope to solve (^9|) or (|3^) by numerical in- 
tegration. However, the implicit dependence of U [x) and 
(3 on cix) via the three-dimensional integral ( ^2[ ) makes 
this a formidable task. Similar considerations hold for u;, 
which has to be computed from the relation L = \u}. To 
simplify things and in particular to reduce the dimension- 
ality of the integrals involved, we pass to spherical coordi- 
nates, X = {x,9,(f)), and expand the Newtonian potential 



* In order to restore a ground state, however, the Thirring 
potential [^, which mimics Newtonian gravitv, has been used 
together with Boltzmann statistics (see e.g. [p5[). 
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in series of real spherical harmonics (see e.g. |28 ) 



1 47r {x\/x'y 

(34) 



with xM x' — minja;, x'} and x A x' — max{x, a;'}. At the 
same time, we formally expand also the relative occupancy 
c: 



C{x) = X! X! ^lrn{x)Yi 



(35) 



i=0 m=-l 



bim{x) is a radial function whose precise form we will have 
to derive. Using the above series, together with the com- 
pleteness relation for our basis set {Yim}, 



(36) 



Yi^{9,(j))Yi,jn,{9,(j)) d cos 9 d(p ^ SwSr, 
one can easily show that 

Uix)=Y,Ulm{x)Yirni9,cj,) 



l.m 



in f [xM x') 



(37) 



uim{x) = j (^^^,y+i fe''n(a; ){x) dx (38) 

Multiplying both sides of (|3^) by Yim and integrating over 
angular variables one obtains for bun the system of integral 
equations 



him [x)^ g{x,9, ^)Yim (9, ct>) d cos 9 dcf) (39) 



g{x,9,cj)) 



where Z = 0, 1, . . . and m = 
uim, P and u) depend on hi 



I, —I + 1, . . . ,1. Notice that 
This system is completely 
equivalent to (|3C|), but at least an iterative solution pro- 
cedure is imaginable. After having fixed E and L, start- 
ing from an initial reasonable guess for bimix), one can 
compute ui„i{x) from ( p8|) (1-dim. integral) (and /? from 
(^)). Using this, 6;„i(xj can be re-calculated from ( ^ 
(2-dim. integral) to improve the guess, and so until con- 
vergence via e.g. a simple Newton-Raphson method. As 
the starting point, it is convenient to take a high-energy 
configuration, where the kinetic contribution is expected 
to be much larger than the gravitational one, and a homo- 
geneous cloud ("gas") type of c{x) is very likely to solve 
(pO|). In particular, one can initiate from a totally sym- 
metric configuration where 600(2;) = and bim{x) ~ for 
(Z,m)^(0,0). 

Clearly, actual calculations must be performed with a 
finite number of harmonics /max) i-e- the series ( p4[ ) must 
be truncated. The effects of such a truncation are shown 
m Fig. |l|. One sees the potential felt by one particle due to 
another one fixed at the position x — 0.4 with = 0. In 
(a), the first particle moves from a; — to a; = 1 at fixed 



Phi=Pl 8 lmax=2 



Phl=Pl 8 lmax=4 





0.2 0.4 0.6 0. 

Distance 



Phi=Pl 8 lmax=8 



0.20.40.60.8 1 

Distance 



Phi=Pi 8 lmax=16 





0.2 0.4 0.6 0. 

Distance 



0.20.40.60.8 1 

Distance 



(a) Fixed <j) = n/S, variable In 



Phi=0 lmax=16 




Phi=Pi 16 lmax=16 



0.2 0.4 

Distance 



Phi=Pi 8 lmax=16 




0.2 0.4 0.6 

Distance 



Phi=Pi 4 lmax=16 




0.2 0.4 0.6 

Distance 



0.2 0.4 0.6 0.8 

Distance 



(b) Fixed imax ~ 16, variable (j) 

Fig. 1. Effects of truncation of the series (^^ to the term of 
order imax (see text for details). Dashed and continuous lines 
represent the truncated potential and the Newtonian potential, 
respectively. 



4> = n/8 and the true potential is compared with the trun- 
cated one, with maximum number of included harmonics 
variable from 2 to 16. The latter case clearly reproduces 
the Newtonian force with a good degree of accuracy. How- 
ever, the 0-dependence must be considered also. This is 
shown in (b), where we fixed Imax = 16 and measured 
the potential varying the (p of the second particle, keeping 
the first one fixed. It is clearly seen that the truncated 
potential works fine sufficiently far from the "probe" par- 
ticle, while small deviations occur when the particles are 
too close. However this short-distance problem is substan- 
tially cured by our choice to deal with non-overlapping 
particles. A further hint about what should be the maxi- 
mum order of harmonics to be included in the calculation 
comes from the study of the behaviour of bimix) for typ- 
ical solutions of (p9|), an example of which is reported in 
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X 

Fig. 2. Typical behaviour of the radial function bim{x) for 
different I {I even, 2 < I < 14) at fixed m = 2. This particular 
plot was obtained for E = —0.18 and L = 0.44. 



Fig. ||. One clearly sees that bim dies out as I increases, 
and that already for Z = 14 it is for all practical purposes 
zero. 

Hence, we solved (|39| ) taking Zmax = 16. We also ex- 
cluded odd harmonics. Simple symmetry considerations 
suggest that exclusion of / = 1 harmonics fixes the center 
of mass in the origin, while absence of higher-order odd 
harmonics prevents the formation of asymmetric struc- 
tures (e.g. two clusters of different sizes lying at different 
distances from the origin). Their effects on the phase dia- 
gram will be studied elsewhere j29j . Finally, we measured 
energy and angular momentum in units of GN'^/R and 
(i?G'iV3)i/2^ respectively, and took O = 0.02 alwayfj^ The 
results of this analysis are reported in the next section. 



3 Thermodynamics 

3.1 Phase diagram 

We shall discuss here solutions of (|3^) obtained at fixed E 
and i in a slab of the {E, L) plane delimited by the lines 
E — L ~ I and E — L. The entr opy corresponding to each 
solution can be calculated via (|27|). Pure thermodynamic 
phases with one (macroscopic) equilibrium state can be 
discerned from phase coexistence regions by studying the 
Hessian of S* in and L, i.e. 

Hes(.,,)[5]^det(^3f/-|f) (40) 

In the microcanonical ensemble pure phases are character- 
ized as having HeS(£; > (the entropy as a function 
of E and L is concave), while phase coexistence regions 

^ The ©-dependence of the results is an important issue. In 
fact, when & is too large even the usual gravitational collapse 
transition does not take place because of particles jamming 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 

E 



Fig. 3. Phase diagram in the {E, L)-plane. The dashed lines 
E — L — 1 (left) and E = L (right) delimit the region where 
the Hessian was calculated. The four markers (x) correspond 
to the four situations described in Fig. H. 



have Hes(£^) [S] < ^,|2^ (the entropy has a convex in- 
truder). Mixed phases are characterized by negative spe- 
cific heat and hence ensemble inequivalence. (The reader is 
referred to for an introductory account on microcanon- 
ical thermostatistics.) Fig. ^ shows the resulting global 
phase diagram of the system. In Fig. ^ one sees a sample 
of stationary distributions, together with the values of E 
and L at which they were obtained. 

Results can be summarized as follows. For low angular 
momenta: 

a. at high energies (the kinetic term dominates) , the solu- 
tion of ( ^9|) is unique and is of the homogeneous cloud 
("gas") type; one finds that HeS(£;^L)[S'] > and the 
corresponding (pure) phase is labeled as "gas" phase; 

b. at low energies, where gravity dominates, the solution 
of (|9|) is unique and of the single-cluster type ( "single 
star", e.g. Fig. ^); correspondingly, one finds a pure 
thermodynamic state with Hes(£;_i) [S*] > 0, which we 
call "single-cluster" state; 

c. in between these two regimes, one finds a phase coexis- 
tence region with ilesi^E,L) [S] < and negative specific 
heat, in which different solutions occur at each point 
{E,L) ("mixed phase"); here, single-cluster and gas 
type of solutions occur. 

For slowly rotating systems one thus recovers a scenario 
that is similar to the usual gravitational collapse that is 
found in theories without angular momentum. For higher 
angular momenta (i.e. for a rapidly rotating system), in- 
stead: 

d. at high energies, the situation gas-type of solutions are 
obtained, and the situation is as in a. above; 

e. at low energies, the Hessian is positive and solutions 
of (|39) are of double-cluster type ("double star", e.g. 
FigT^b), and the corresponding phase, labeled "double 
cluster" is pure; 
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(c) E = -0.06, L = 0.4 (d) E = -0.42, L = 0.5 



Fig. 4. Examples of stationary distributions c[x) occurring 
inside our spherical box. Shown are the contour plot and, above 
it, the density profile: (a) "single star", (b) "double star", (c) 
"disk", (d) "ring". 



f. at intermediate energies, multiple solutions are found, 
both of double-cluster type and deformed gas-type. 
The latter in particular can be disks (e.g. Fig. ^) or, 
if the angular momentum is high enough, rings (e.g. 
Fig. ^). Correspondingly, the Hessian is negative and 
we have a phase coexistence region with negative spe- 
cific heat. 

The system thus turns out to have three pure phases 
("gas", "single star" and "double star"), separated by a 
large mixed phase. The occurrence of double-star-like so- 
lutions is the most remarkable effect of rotation. To get an 
idea of the coexistence of different solutions in the mixed 
phase, in Fig. ^ we plot the entropy for three different so- 
lutions (ring, single cluster, double cluster) in a range of 
energy at fixed L = 0.5. One sees that the entropy has 
a "convex intruder" , corresponding to the mixed phase 
and implying negative specific heat. The three solutions 
coexist in a whole range of energies, while at low ener- 
gies, as evident from the phase diagram, the double-cluster 
solution only survives. The point where the rotationally 
asymmetric double-cluster solution bifurcates from the ro- 
tationally symmetric ring solution corresponds to the be- 
ginning of the mixed phase at L = 0.5. In the mixed phase: 
double-cluster and (deformed) gas type of configurations 



2.4 I • • • • ' • • • • ' • • • • 1 

-1 -0.5 0.5 

E 

Fig. 5. Entropy as a function of energy at fixed L = 0.5 for 
three different solutions, as shown. To make the convex in- 
truder in the entropy (corresponding to the mixed phase) more 
evident, we subtracted the quantity (15/4)i5 from the entropy. 



compete at high E and L; single-cluster and gas compete 
at low L; finally, single-cluster and double-cluster compete 
at intermediate values of L. 

The crucial issue of stability (i.e. which of these config- 
urations are actually entropy maxima in the space of c's) 
will be dealt with in Sec. 3.4. For the moment, let it suffice 
to say that in the mixed phase, rotationally symmetric 
structures are unstable to perturbations that break ro- 
tational symmetry. Hence ring configurations, which also 
occur in the mixed phase, are not stable. Deformed gas 
configurations (e.g. disks or rings) that occur in the "gas" 
phase, e.g. close to the phase boundary, are stable. Be- 
fore analyzing the different transitions that take place, we 
shall briefly discuss the important issue of negative specific 
heat. 



3.2 Caloric curves, specific heat 

Once the solutions are obtained, the microcanonical en- 
tropy surface can be immediately calculated from ( p7| ) . It 
is reported in Fig. ^ together with the /3-surface, namely 

P^PiE,L)=(^) ^1 (41) 

\ / L=constant ^ 

representing the inverse microcanonical temperature as a 
function of energy and angular momentum. The central 
region in the entropy surface corresponds to the mixed 
phase and has HeS(£;^i)[S'] < 0. Slices of the /3-surface at 
different L (caloric curves) are shown in Fig. ^. One sees 
that [3 is increasing with E for a whole range of energies. 
This means that in that range ^ > 0, or equivalently 
that ^ < 0, i.e. that the specific heat is negative. Physi- 
cally, if the system is heated (increase of total energy) its 
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4.5 



(b) 

Fig. 6. (a) Microcanonical entropy 5* as a function of E and 
L. In the central region the Hessian Hes(B_L) [S] is negative, 
signaling the separation of multiple phases (mixed phase). In 
the three clear regions, the Hessian is positive and pure ther- 
modynamic phases (gas, single star, binary star, respectively) 
occur, (b) Behaviour of /3 = OeS as a function of E and L. 



-0.4 -0.2 0.2 0.4 0.6 

E 



(c) L = 0.6 

temperature diminishes, and vice-versa if it loses energy 
its temperature increases. This is well-known to happen in 

stars. As they irradiate and release energy, they become ^'S- ^''^^^ /3-surface (inverse microcanonical 

1,, jiij. J iioi,!.! - V. temperature, caloric curves) at diuerent angular momenta, 

hotter and hotter and contract. Such a behaviour, how- t- ^ ; & 

ever, is not peculiar to self-gravitating systems, but is the 
generic signal of a phase separation in the microcanonical 
ensemble of finite systems and is not recoverable in 

the canonical ensemble, where the specific heat is propor- 
tional to energy fluctuations, i.e. non-negative definite. 
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3.3 Phase transitions, symmetry breaking 

We now turn to investigating the phase transitions. The 
low-angular-momentum collapse transitions are analogous 
to those discussed at length in the literature. A probably 
convenient order parameter to describe them is the density 
contrast, defined as the center-to-edge ratio of the parti- 
cle density. We shall concentrate here on the rotational- 
symmetry-breaking transition to double-cluster solutions, 
which is the truly novel phenomenon introduced by rota- 
tion. A convenient order parameter to detect such a tran- 
sition is 



(42) 




-0.8 -0.6 -0.4 -0.2 



that is, the difference between 1 and 2 diagonal compo- 
nents of the inertia tensor (|o|). One expects D to be 
zero when the solution is rotationally-symmetric, and non- 
zero for a solution without rotational symmetry. The rea- 
son is physically clear. If L = (i.e. in absence of ro- 
tation) the system is necessarily isotropic {In = I22 — 
and rotational symmetry cannot be broken. When 
L ^ 0, anisotropics may occur (133 ^ Iii,l22) and one 
can have either rotationally-homogeneous (Jn = 122) or 
rotationally-heterogeneous (/n 7^ 122) solutions. The lat- 
ter correspond to double clusters. (We remind the reader 
that the angular momentum is chosen to lie parallel to 
the 3-axis). In Fig. || we show explicitly that D actu- 
ally behaves as a conventional order parameter for the 
rotational-symmetry-breaking transition (and the appear- 
ance of double clusters) by plotting it as a function of E at 
low and high angular momentum. By comparing Fig. ||b 
with the phase diagram the reader can notice that the 
transition occurs exactly at the phase boundary. 

Another view of the same transition is given in Fig. ^. 
One sees the microcanonical entropy as a function oi In — 
I22 at fixed L = 0.5 and different energies. The four di- 
agrams shown correspond to the four markers displayed 
across the gas-double cluster phase boundary in the phase 
diagram. The entropy is clearly seen to develop two peaks 
at non-zero values of In — J22, corresponding to double 
clusters systems, with the two stars either aligned on the 
1-axis or on the 2-axis, respectively. The fact that S be- 
comes fiat at the phase transition indicates that the tran- 
sition is second-order. The minimum of S, occurring at 
III = -^22 7 corresponds to another, rotationally-symmetric 
solution of (|3C|). In particular, it is a ring. This brings us 
to the problem of stability and clarifies further the struc- 
ture of the mixed phase: at fixed E and L, the entropy in 
the c space has (at least) two maxima corresponding to 
double stars aligned on different axes. These are the only 
stable configurations. 



(a) L = 0.2 



(b) L = 0.6 

Fig. 8. Behaviour of the order parameter D = Jn — I22 as a 
function of energy E at fixed angular momentum L. At high 
L, breaking of rotational symmetry is signaled by a D 7^ 0. 



the reference frame of the principal inertia axes, where I 
is diagonal, one can calculate such a variation explicitly. 
Omitting details of the lengthy calculation, one finds 



13 f Sc{x)Sc{x') 



0^ J \x — X 
2 

^ 302 



dxdx'— — 



{5c{x)f 



5c{x) log T—--^^dx 



ej c{x){i-c{x)) 

2 



dx- 



1 - c{x) 



3 

Jaa \_J 



Sc{x] 



(43) 



3.4 Stability 

Usually, the analysis of the (local) stability properties of 
the stationary points of the microcanonical entropy ( p7| ) 
at fixed mass, energy and angular momentum relies on 
the study the sign of second variation of the entropy. In 



with Sc{x) a mass-preserving perturbation (J 5c(x)dx = 
0). I*^"^ = l'^°^[c] stands for the a-th column of the inertia 
tensor L It is easy to show that for (3^0, that is at 
sufficiently high energy where the kinetic term dominates, 
homogeneous gas-like stationary configurations are stable 
against any perturbation. It suffices to observe that when 
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-0.2 -0.1 0.0 0.1 -0.4 -0.2 0.0 0.2 0.4 



2.644 



2.641 - 



2.638 
2.460 



2.457 - 



2.454 







E=-0.04 






£=-0.7 



2.555 

2.550 

2.545 

2.540 
2.37 

2.36 

2.35 



-0.4 -0.2 0.0 0.2 0.4 -0.4-0.2 0.0 0.2 0.4 

Fig. 9. Entropy as a function of /n — I22 at L = 0.5 and 
different values of E. The values of E and L for the four figures 
correspond to the four markers (x) shown in Fig. ^ 



/3 ^ the above equation reduces to (k > constant) 

<5^S'(/3 ^ 0) = -K / {5c{x)fdx < (44) 



because c{x) — constant and 6c(x) is mass-preserving by 
assumption. 

In a more general setting, the second variation of the 
entropy can be written as the quadratic form 



S^S^ J f{x)K{x,x')f{x') dx dx' 

where the kernel K is given by 

(52<? 1 5{x - x') 



(45) 



K{x,x') = 

2/3^ r 5<P 
3 [5c(a;) 

U} 

2 



5c{x)5c{x') c{x){l - c{x') 



1 rr. SI 

LJ UJ 

2 ' ' ' 
SI 



-I 



Sc{x) 

SI 



Sc{x') 

51 , 



1 rj. SI 

r 



SI 



Sc{x') 



Sc{x') Sc{x) Sc{x) Sc{x') 



IV (46) 



and where as before <P stands for the Newtonian poten- 
tial. Stability analysis is then equivalent to studying the 
eigenvalue problem for K, namely 



K{x,x')f{x')dx' = \f{x) 



(see e.g. 30|). From a mathematical viewpoint, this task is 
extremely sophisticated. We shall therefore limit ourselves 
here to discuss the stability of rotationally-symmetric con- 
figurations against perturbations that break rotational sym- 
metry, deferring the reader to a later publication for a 
more complete analysis of the stability problem. 



From the analysis of the preceding section, and in par- 
ticular from Fig. ^, it stems that rotationally symmet- 
ric configurations become unstable against perturbations 
that break rotational symmetry at sufficiently high an- 
gular momenta and energies. Two rotationally asymmet- 
ric solutions of ( p9| ) bifurcate continuously from the rota- 
tionally symmetric state. At least for what concerns this 
class of perturbations, it is safe to claim that rotationally 
symmetric structures are stable up to the phase bound- 
ary. Among them, one finds gas-like homogeneous config- 
urations, and deformed-gas configurations (i.e. disks and 
rings). At the phase boundary, solutions with D — be- 
come entropy minima at least along one direction in the 
c space, and are no longer stable against rotational sym- 
metry breaking. Instead, the only stable configurations in 
the mixed phase at high enough angular momentum are 
double-cluster like {D ^ 0). 

A more technical argument that supports this conclu- 
sion is the following. Let (see (|30|)) 



G[c] = (l + ( 
It is easy to show 
S^S 



,#t/(x)-i/3(u,xx)^-t-/.^l-l 



)--c{x) 



that 



Sc{x)Sc{x') 



1 



Gfcl=0 



SG 
Sc{x') 



(48) 



(49) 



G[c]=0 



with 7 > 0, i.e. that the second functional derivative of 
the entropy evaluated at the stationary point vanishes to- 
gether with the functional derivative of G evaluated at 
the same point, and that the two have the same sign. 
This implies that the stability analysis can be reduced 
to the study of the sign of SG/Sc. However, the latter 
problem is dealt with when applying the Newton- Raphson 
method to solve (|3C|). Starting from high energies, in order 
to provoke a rotationally asymmetric solution an appro- 
priate external field (in this case, it is related to the order 
parameter discussed in the previous subsection) must be 
added to G as an external perturbation. A bifurcation 
of a rotationally-asymmetric solution implies a change of 
sign of SG/Sc, meaning that the rotationally symmetric 
one has become unstable to that particular perturbation. 
Hence, the instability-onset line for perturbations that 
break rotational symmetry numerically coincides with the 
phase boundary between the "gas" and the "double clus- 
ter" phases, where Hes^^; — 0, displayed in the phase 
diagram. 



4 Outlook and final remarks 



We have presented an analysis of the equilibrium prop- 

(47) 

erties of a self-gravitating and rotating gas using a mi- 



It is sufficient to note that 
SS 



Sc{x) 



.llogll 

6> ^ c(a;)(l 



c(x))(c(x) + G[c]) 



cix)-G[c]) 



Taking the functional derivative of this expression with respect 
to c{x') and evaluating it at the stationary point, one finds ( |49| ) 
with 7 = [0c{x){l - c{x))]-K 
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crocanonical mean-field approach. Our main result con- 
cerns the spontaneous breaking of the rotational symme- 
try, which takes place at high angular momentum and 
gives rise to non-trivial density profiles, e.g. "double stars" : 
a rapidly rotating iV-body system kept together by gravi- 
tation only at equilibrium spontaneously organizes in two 
distinct dense clusters, provided the ratio between rota- 
tional and gravitational energy is sufficiently high and the 
total energy is not too large. We have derived the global 
phase diagram of the model and discussed the related 
thermodynamic picture, providing a phenomenological de- 
scription of the phase transitions occurring and analyzing 
the stability of high-energy rotationally-symmetric equi- 
librium states against perturbations that break rotational 
symmetry, showing that non-trivial rotationally symmet- 
ric solutions such as rings become unstable in the mixed 
(phase coexistence) region. To our knowledge, these re- 
sults constitute the most complete equilibrium description 
of a self-gravitating and rotating system to date. To con- 
clude, we would like to put forward some final remarks 
and open problems. 

i. While it is certainly possible to improve on the results 
presented here by increasing the maximum order of 
the even harmonics included in the calculation, we do 
not expect any major qualitative difference with the 
picture we describe here (Zmax — 16). The inclusion of 
odd harmonics would instead lead to the formation of 
asymmetric double clusters, and to a more complete 
phase diagram and a full classification of the different 
possible equilibrium configurations as a consequence. 
This issue will be treated elsewhere [ p9[ . 

ii. Our results stem from a mean-field analysis, in which 
particle-particle correlations are completely neglected. 

iii. We have not studied how our results depend on 0, 
namely on the average density of the system. We men- 
tioned that this is an important issue, certainly worth 
to be investigated with great care. 

iv. A general stability analysis requires, as we said, the 
study of the eigenvalue problem for K, Eq. (^). The 
most interesting open question is that of marginal sta- 
bility, that is solutions of ( p7| ) with zero eigenvalue. 
However, the approach we discussed in the last sec- 
tion, connecting the stability analysis to the bifurca- 
tion analysis, describes correctly the stability of the 
different states against perturbations that break rota- 
tional symmetry (within the included number of har- 
monics) . 

V. From a physical point of view, it is known that en- 
ergy and angular momentum are possibly not the only 
conserved quantities to be taken into account if one 
wants to recover some observational features of e.g. 
galaxies | |3T| . 

vi. Of course, we have dealt with equilibrium properties 
exclusively, and provided a kind of classification of the 
different possible states of the system. This clearly 
leaves open many important questions concerning dy- 
namics and relaxation mechanisms, which are believed 
to be particularly subtle in self-gravitating systems 

0- 



On a more general level, the equilibrium properties of 
non-extensive Hamiltonian systems (self-gravitating and 
rotating systems being the most important example) can 
be well described by Boltzmann's principle (^) . We would 
finally like to stress the potentialities of the method pre- 
sented here for studying the equilibrium properties of sys- 
tems with long-range (see footnote 1) forces. Its basic in- 
gredients are: (a) microcanonical statistics and (b) a series 
expansion of the potential using a proper basis set. This 
same scheme can easily be modified to work with other 
types of long-range potentials (e.g. Coulomb or Yukawa). 
In the light of this last remark, self-gravitating and rotat- 
ing systems (for which the use of "corrective" techniques 
such as the Kac-Uhlenbeck-Hemmer method js^] is sub- 
stantially not helpful in clarifying its complex thermody- 
namic structure) should be considered as a fundamental 
testing ground for techniques to analyze the statistical 
equilibrium properties of systems with long-range inter- 
actions. 
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